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ABSTRACT 


A  linear  stability  analysis  of  the  Benard  problem  for  two  layers  of  different  fluids  lying 
on  top  of  each  other  and  bounded  by  free  surfaces  is  considered.  The  fluids  are  assumed  to 
be  similar  and  perturbation  methods  are  used  to  calculate  the  eigenvalue  in  closed  form. 
The  case  of  the  Rayleigh  number  and  wavenumber  of  the  disturbance  being  close  to  the 
first  criticality  of  the  one-fluid  Benard  problem  has  been  investigated  in  a  previous  paper  , 
and  was  found  to  exhibit  both  overstability  and  convective  instability.  In  this  paper,  the 
Rayleigh  number  is  assumed  to  be  less  than  that  of  the  first  criticality  of  the  one-fluid 
problem,  and  in  this  situation,  overstability  does  not  occur.  An  unexpected  result  is  that 
by  an  appropriate  choice  of  parameters,  it  is  possible  to  find  linearly  stable  arrangements 
with  the  more  dense  fluid  on  top. 
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T: 


SIGNIFICANCE  AND  EXPLANATION 

Interfacial  stability  of  a  two-layer  convective  problem  is  analyzed  for  the  case  of  two 
similar  liquids.  Applications  may  involve  the  modeling  of  convective  instability  in  a  variety 
of  multi-layered  films,  for  example,  pattern  formation  in  rocks. 
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INTERFACIAL  STABILITY  IN  A 
TWO-LAYER  BENARD  PROBLEM 

Yuriko  Renardy 


I.  INTRODUCTION 

Two  layers  of  fluids  are  sandwiched  between  infinite  parallel  horizontal  boundaries. 
The  fluids  have  only  slightly  different  properties.  We  neglect  any  diffusion  of  one  fluid  into 
the  other,  so  that  there  is  an  interface,  with  surface  tension  taken  into  account.  The  lower 
plate  is  kept  at  a  slightly  higher  temperature  than  the  upper  plate,  and  the  tangential 
stress  and  normal  velocity  vanish  at  the  boundaries.  Although  the  “free-free”  boundary 
conditions  are  physically  unrealistic,  they  yield  the  advantage  that  closed-form  solutions 
can  be  obtained,  and  their  analysis  helps  to  shed  light  into  the  numerical  investigation  of 
the  more  realistic  “rigid-rigid” boundary  problem2.  We  consider  the  linear  stability  of  the 
rest  state  when  the  Rayleigh  number  is  below  the  first  criticality  of  the  corresponding  one- 
fluid  problem.  The  case  when  the  Rayleigh  number  and  wavenumber  of  the  disturbance 
are  close  to  that  of  the  first  criticality  of  the  one-fluid  problem  has  been  investigated 
previously1.  For  the  linear  stability  analysis,  we  only  need  to  consider  two-dimensional 
disturbances,  because  of  the  rotational  symmetry  about  the  vertical  axis.  The  questions 
we  ask  are:  what  kinds  of  instabilities  can  arise;  are  they  different  or  similar  to  those  that 
arise  in  the  one-fluid  problem;  and  how  do  these  instabilities  depend  on  the  dimensionless 
parameters. 

The  equations  used  in  each  fluid  are  the  Navier-Stokes  equations  with  the  Boussinesq 
approximation,  the  linear  heat  equation  and  incompressibility.  Temperature  differences 
in  the  problem  are  assumed  to  be  small.  The  properties  of  the  fluids  are  assumed  to  be 
constants  except  for  the  density  in  the  gravity  term  in  the  Navier-Stokes  equations,  where 
the  density  is  expanded  as  a  Taylor  series  about  the  temperature  of  the  top  boundary,  and 
truncated  so  that  it  is  a  linear  function  of  the  temperature. 


At  the  interface,  the  following  conditions  are  assumed  to  hold:  the  kinematic  free- 
surface  condition,  the  continuity  of  velocity  and  shear  stress,  the  difference  in  the  normal 
stresses  must  be  balanced  by  surface  tension,  and  the  temperature  and  heat  flux  are 
continuous.  These  equations  yield  9  dimensionless  parameters:  a  Rayleigh  number  and 
Prandtl  number  based  on  one  of  the  fluids,  a  surface  tension  parameter,  and  6  ratios  of 
the  various  thermal  and  mechanical  properties  of  the  fluids. 

One  solution  to  the  equations  is  the  rest  state,  with  a  linear  temperature  gradient  in 
each  fluid,  and  a  flat  interface.  We  consider  the  linear  stability  of  this  solution  by  adding 
a  small  perturbation  which  is  proportional  to  exp(iax  +  at)  where  x  is  the  dimensionless 
horizontal  variable  and  t  is  the  dimensionless  time.  The  problem  is  set  up  as  an  eigenvalue 
problem,  to  calculate  o  in  terms  of  the  other  parameters. 

It  is  instructive  to  recall  some  of  the  results  for  the  one-fluid  problem.  Rayleigh, 
in  1916, 3  solved  the  linear  stability  problem  for  one  fluid,  for  “free-free”  boundaries.  The 
eigenvalues  and  eigenfunctions  were  found  in  closed  form,  and  the  marginal  stability  curves 
are  given  by  Rj  -  (j27r2  +  o2)3/a2,  j  =  1,2,....  R\  <  R.2  <  Rz  <  ...  The  critical  Rayleigh 
number  is  27n*/4.  The  Prandtl  number  does  not  enter  into  the  criticality  condition,  but 
appears  in  the  growth  rate.  The  “rigid-rigid”  boundary  and  “free-rigid”  boundary  cases 
were  solved  later3:  here,  the  eigenvalues  and  eigenfunctions  cannot  be  found  in  closed  form 
but  have  to  be  calculated  numerically. 

Another  well-known  result  for  the  one-fluid  problem  is  the  “exchange  of  stabilities”: 
the  equations  are  self-adjoint  so  that  all  the  eigenvalues  are  real.  A  consequence  is  that 
the  velocity  at  marginal  instability  is  not  periodic.  The  exponential  growth  of  small  dis¬ 
turbances  leads  to  a  steady  nonlinear  motion,  for  example,  steady  cellular  convection. 
However,  in  the  presence  of  solutes  with  one  or  more  concentration  gradients,  or  continu¬ 
ous  density-stratification,  a  time-periodic  marginal  instability  (  “over stability”)  is  possible. 

In  the  two-fluid  problem,  the  equations  are  not  self-adjoint2.  Therefore,  overstability 
is  possible,  and  such  a  situation  was  found  numerically  for  the  case  of  rigid  boundaries. 
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The  marginal  eigenvalues  are  a  complex  conjugate  pair  of  multiplicity  two.  The  rest  state 
would  bifurcate  to  either  a  traveling  wave  or  a  standing  wave  solution.  An  analysis  of  the 
subsequent  nonlinear  problem  may  be  related  to  that4  of  one-fluid  double  diffusion. 

It  would  be  interesting  to  find  out  under  what  conditions  complex  eigenvalues  arise 
and  what  conditions  yield  real  eigenvalues  in  the  linear  stability  problem,  but  a  numerical 
attempt  at  this  is  not  feasible  because  of  the  large  number  of  parameters.  However, 
the  case  of  free  boundaries  where  the  two  fluids  have  thermal  and  mechanical  properties 
that  differ  by  a  small  amount  of  0(c)  can  be  analyzed  with  perturbation  methods  in 
the  parameter  c.  The  leading  terms  in  the  perturbation  expansion  for  the  eigenvalue  are 
found  in  closed  form.  In  a  previous  paper1,  this  analysis  was  performed  for  the  situation 
where  the  unperturbed  one-fluid  problem  is  close  to  the  first  criticality.  Thus,  the  analysis 
involved  the  perturbation  of  a  double  zero  eigenvalue. 

In  the  present  paper,  the  unperturbed  one-fluid  problem  is  below  the  first  criticality,  so 
that  the  analysis  for  two  similar  fluids  involves  the  perturbation  of  a  simple  zero  eigenvalue. 
Since  the  equations  are  real,  the  perturbed  eigenvalue  is  also  real.  Hence,  when  the  interface 
is  unstable,  one  would  expect  the  fluids  to  go  into  a  different  arrangement.  For  example, 
if  the  instability  results  from  the  upper  fluid  being  the  heavier,  one  would  expect  that 
it  should  fall  and  the  fluids  should  exchange  positions.  However,  there  are  more  subtle 
cases  as  will  be  discussed.  The  short-wave  asymptotics  for  a  is  identical  to  that  obtained 
for  the  rigid-rigid  boundaries2,  and  indeed  this  agrees  with  the  short-wave  asymptotics 
of  our  perturbation  formula.  Differences  in  density,  coefficients  of  cubical  expansion  and 
surface  tension  are  important  in  the  short-wave  limit.  In  the  long-wave  limit,  the  volume 
ratio,  and  differences  in  thermal  conductivity,  density  and  coefficients  of  cubical  expansion 
are  important.  By  making  appropriate  choices  for  the  parameters,  it  is  possible  to  find 
a  linearly  stable  arrangement  with  the  more  dense  fluid  on  top.  In  fact,  the  formula 
for  a  for  small  /j  shows  the  unexpected  “thin-layer  effect",  where  stability  is  strongly 
influenced  by  the  thermal  conductivity  stratification.  This  behavior  is  reminiscent  of  the 
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lubricating  effect  of  viscosity  stratification  in  parallel  shear  flows  composed  of  two  layers 
of  different  fluids5,  where  linearly  stable  arrangements  with  the  heavier  fluid  lying  on  top 
have  been  found,  provided  that  suitable  values  are  chosen  for  the  volume  ratio,  surface 
tension,  viscosity  difference  and  density  difference. 
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II.  GOVERNING  EQUATIONS 

We  consider  two  fluids,  lying  on  top  of  each  other  between  two  parallel  boundaries 
of  infinite  extent  in  the  (x‘ ,  z")-plane.  Fluid  1  denotes  the  lower  of  the  fluids  and  Fluid 
2  denotes  the  upper  fluid.  Asterisks  denote  dimensional  variables.  The  upper  boundary 
at  r'  =  /‘  is  kept  at  a  constant  temperature  Tt-,,  and  the  lower  plate  is  kept  at  a  higher 
constant  temperature  T(,  -f  A  T'.  The  interface  is  at  z  =  l\.  Subscript  i  (=1,2)  denotes 
Fluid  i.  At  temperature  T^.  the  fluids  have  coefficients  of  cubical  expansion  d,  ,  thermal 
diffusivity  #c,,  thermal  conductivity  fc,,  viscosity  pt ,  kinematic  viscosity  j/,  and  density  pt. 
We  use  the  Boussinesq  approximation  in  the  Navier-Stokes  equations,  that  is,  the  densities 
in  the  buoyancy  term  for  each  fluid  is  approximated  by  p,(l  -  d,(T'  -  T0')). 

There  are  6  dimensionless  ratios  of  the  physical  properties: 


m  = 

r  = 

1  = 

?  = 


M2’ 

Pi 

> 

P  2 

*2 

h. 

h' 


(i) 


a, 

a2 

/.  =  =  1  -  h- 

Following  Drazin  and  Reid  3,  we  introduce  dimensionless  variables  (without  asterisks) 
as  follows: 

(x,z)  =  (x,z')/r, 

t  =  Kit'/r2,  (2) 

U  =  uT  j K J, 

T  =  T'/AT’, 
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p=  p'/'2/(pi«?)- 

Here,  u‘  =  (u’,u;*)  is  the  velocity,  p’  the  pressure  and  T‘  the  temperature.  We  define 
a  Rayleigh  number 

R  =  gaiAT'l^/iKiVt), 

where  g  denotes  gravitational  acceleration,  and  a  Prandtl  number 

P  =  U\  j  K\ . 

We  include  a  surface  tension  between  the  fluids,  described  by  a  dimensionless  parameter 
S  =  5*/"/(acj/zi),  where  S’  is  the  dimensional  surface  tension  coefficient. 

We  study  the  linear  stability  of  the  static  solution 

u  =  0, 

T  =  T0+1-  Aiz  for  0  <  2  < /,,  (3) 

=  To  +  -42(1  ~  2)  for  li  <  z  <  1, 
where  Ai  =  and  A2  =  gAi. 

If  disturbances  are  proportional  to  exp(cr/  +  tax),  then  the  following  linearized  eigenvalue 
problem  arises  for  the  velocity  u  and  perturbation  0  to  the  temperature. 

For  0  <  2  <  /  j , 


oQ  =  wA  1  +  A0, 
ou  =  -dp/dx  +  P  A  u. 
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aw  =  -dpfdz  +  RPQ  +  PA11, 


du  dw 
dx  dz 


0. 


(4) 


For  li  <  z  <  1, 


(70  =  wA  2 


-  A0. 

1 
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dp  r  _  . 

ov  =  —  r  — — l-  — P  A  u, 

<9r  m 

dp  RP  _  r  . 

era’  =  -r  — ^ -  +  — —  ©  -1-  —  P  A  w,  (5) 

dz  p  m 

du  dw 
¥x  ~  dl  ~  ' 

At  z  =  0,1,  slip  boundary  conditions  are  assumed  to  apply,  i.e.,  the  shear  stress  is 


zero,  or,  =  0: 


_  du 

e  =  w  =  —  =  o. 

dz 


The  interface  is  perturbed  to  the  position  z  =  l\  +  h(x,t).  The  interface  conditions 
linearized  about  z=/j  are2,G  : 

continuity  of  temperature:  J©  J=  h  JA  |, 
continuity  of  heat  flux:  j=0, 

continuity  of  velocity:  flw  J=  [u  J=0,  (7) 

continuity  of  shear  stress:  |/z(|^  +  §7)  J=0, 
balance  of  normal  stress: 

,dw\  1  dw  9.  d2h 

p,-pl+2P^---^)  +  Mlk-rS^=0, 

kinematic  free  surface  condition: 

oh  —  W\.  (8) 

Here,  j.  Jdenotes  the  jump  of  a  quantity  across  the  interface,  for  example. 


=  Ai  -  A 2, 


and  we  have  set 


.*  V  / 


We  now  proceed  as  follows.  It  is  known2,3  that  when  the  fluids  have  identical  prop- 
rties,  the  eigenvalues  are  real,  the  presence  of  the  interface  introduces  a  zero  eigenvalue 
t  any  Rayleigh  number,  and  another  zero  eigenvalue  occurs  first  at  R  =  ~n*  (or  a  = 
"  *'  27r.  We  let  R  be  arbitrary  but  less  than  the  value  of  the  first  criticality  of  the  one-fluid 
roblem.  We  assume  the  fluids  have  slightly  different  properties,  and  introduce  a  small 


arameter  c.  We  regard  1  -  m,  1 

-  r,  1 

1 

l 

T* 

1 

0,Mi 

and  5  as  small  quantities 

roportional  to  c;  that  is,  we  set 

1 

-  m  =  tfu. 

1 

-  r  =  re. 

1 

~  1  = 

1 

1 

II 

1 

Mi  —  M\t, 

Mx  = 

+  h[r 

5  =  St. 

At  c  =  0.  there  is  a  simple  eigenvalue  a  —  0,  arising  from  the  presence  of  the  interface. 
3r  small  t,  this  eigenvalue  is  perturbed  into  a  power  series  in  f.  The  purpose  of  the 
Mowing  analysis  is  to  find  the  coefficient  of  t  in  this  expansion. 


III.  PERTURBATION  ANALYSIS 


The  perturbation  expansion  for  a  double  eigenvalue  is  carried  out  in  detail  in  Ref  1, 
and  closely  parallels  the  procedure  adopted  here  so  that  the  explanation  in  this  section 
will  be  terse.  The  perturbation  expansion  for  the  simple  eigenvalue  involves  finding  the 
eigenspace  belonging  to  the  eigenvalue  a  =  0  for  the  unperturbed  problem  and  its  adjoint, 
and  does  not  require  finding  the  eigenspace  of  the  perturbed  O(t)  problem  at  all. 

Suppose  o0  is  a  simple  eigenvalue  of  a  matrix  L((.  Let  A  be  an  eigenvector  of  L0 
with  eigenvalue  oo,  and  let  C  be  an  eigenvector  of  L[,  (the  adjoint  of  Lo)  with  eigenvalue 
do  (the  overbar  here  denotes  the  complex  conjugate).  Let  Lo  be  perturbed  into  L{e)  = 
Lo  +  (L\  +  0(f2)  with  L\  depending  smoothly  on  c.  Then  the  peiturbed  eigenvalue  o 
is  given  by  the  zero  of  the  expression  which  represents  to  O(c)  the  projection  of 

L(()  —  o,  first  onto  the  eigenspace  of  the  unperturbed  problem  and  then  onto  the  adjoint 
eigenspace: 

*M)  =  (C,  (Lo  +  (Li  -  o)A)  +  0(c2).  (9) 

For  the  same  reasons  as  in  §111,  Ref.  1,  we  will  need  to  redefine  ty. 

Let  X  denote  the  set  of  functions  (G,u,w,h).  We  introduce  an  inner  product  by 

/•2ir/a  rli 

(xl.x2)  =  /  e,e2  +  U)U2  +  W\W2  dzdx 

Jo  Jz= 0 

r 2*/a  r 1 

+  /  /  9j©2  ■+  UjUo  T  U>lU>2  dzdx 

Jo  J  z  —  li 

r2ir/a  _ 

+  /  hih2  dx  (10) 

Jo 

to  generate  a  Hilbert  space.  In  this  Hilbert  space,  we  consider  the  subspace  determined  by 
the  “Hodge  projection”  (see  space  H  in  Theorem  1.4,  Ref.  7),  that  is,  by  the  conditions  that 
the  velocity  field  be  divergence-free,  that  the  vertical  velocity  vanish  at  the  boundaries,  and 
be  continuous  across  the  interface.  By  L(<)X  we  denote  the  right  hand  sides  of  equations 
(4), (5)  and  (8).  We  regard  L(f )  as  an  operator  in  the  subspace  so  that  the  conditions  on  w 


I  k  *•>’*.  X.' \  i-  *. 


.  r.  v  x*  i.*v  i  *.  ' 


vVOh.’  V  v  ^  v -^  •  -- 


in  (6)  and  (7)  and  the  normal  stress  balance  in  (7)  are  an  integral  part  of  the  definition  of 
L(c).  The  domain  of  definition  of  L(c)  is  determined  by  the  rest  of  the  boundary  conditions 
in  (6)  and  (7),  which  we  write  in  the  form  B(c)X=0.  The  range  of  the  operator  L(c)  must 
satisfy  the  following  conditions  in  order  for  the  pressure  p  occuring  on  the  right  sides  of 
(4)  and  (5)  to  be  determined  as  a  function  of  X:  The  “velocity  part”  of  L(e)X  must  be 
divergence  free,  the  vertical  velocity  must  vanish  on  the  walls  and  be  continuous  across 
the  interface,  and  the  jump  in  p  across  the  interface  must  be  given  by  the  normal  stress 
balance.  Thus,  the  problem  we  wish  to  solve  is:  for  small  c,  find  o  satisfying 


L{e)X  =  oX, 


(11) 


B(e)X  =  0, 

L(c)  =  Lq  +  cLj  +  0(t2), 
B(c)  =  Bo  +  eBi  +0(c2). 


The  explicit  forms  of  L0,  Lx,  B0  and  B\  are  given  in  §111,  Ref.  1. 
We  now  redefine  91: 


*M)  =  (C,  ((L(c)  -  I)"1  -  d)A)  +  0(f2), 


(12) 


where  the  C  and  A  are  as  before. 

The  boundary  value  problem  adjoint  to  (11)  is  calculated  in  Appendix  A  of  Ref.  1. 
Then  we  determine  the  eigenvectors  A  and  C  in  the  Appendix  of  this  paper.  They  satisfy: 


LqA  =  0,  BqA  —  0, 


L0C  =  0,  B'0C  =  0. 


(13) 


In  formula  (12),  we  must  determine  the  expressions 


<C,(L(c)-1)-M) 


(14) 
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to  first  order  in  i.  To  facilitate  this  calculation,  we  introduce  x(l  and  x1  defined  by 


{L{e)-l)-lA  =  x°  +  cxl  +  0(<2).  (15) 

Equating  the  coefficients  of  equal  powers  of  £,  we  find  the  equations  governing  x°  and  x1 

(L0  -  l)x"  =  A. 

B0x 0  =  0,  (16) 

and 

Lji0  +  (Lc,  -  l)x*  =  0, 

B1x°  +  B0x1  =  0.  (17) 

From  (16),  we  find  x°  =  —A.  We  will  not  need  the  solutions  x1  to  the  perturbation  problem 
(17)  but  only  certain  inner  products  involving  them,  namely  (C,x!).  This  is  seen  from 
(12)  and  (15): 

=  (C,x°)  +  e(C,xl)  -  d(C,  A)  +  0(e2).  (18) 

We  calculate  (C,!1)  from  (17): 


(C,x‘)  =  <C,L,x°>  +  (C.Lox1), 

and  an  integration  by  parts: 

(C,  Lox1)  =  (L'0C,x')  +  T, 


(19) 


(20) 


where  the  boundary  integrals  T  are  evaluated  using  the  second  part  of  (17).  (The  boundary 
integrals  would  vanish  if  B ox1  were  zero.) 

We  note  that  x°  =  —A  and  since  LoA  =  0,  =  |^'  =  0  so  that 


L\A 


l  0  \ 
_  M 

dz 

-Vi 

0 
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IV.  RESULTS  AND  DISCUSSION 


We  note  that  a  is  real.  It  is  independent  of  the  Prandtl  number  of  the  basic  one-fluid 
flow,  the  perturbation  on  the  viscosity,  and  the  perturbation  on  thermal  diffusivity.  The 
coefficient  of  f  in  o  is  antisymmetric  about  /j  =  0.5  (see  Figure  l);  the  coefficients  of  M\/P 
and  5  are  symmetric  about  / j  =  0.5  (see  Figures  2-4).  In  M\  IP.  dj  AT'  is  small  compared 
with  1.  In  numerical  computations,  we  have  set  diAT'  =  0.001.  Hence,  the  coefficient  of 
f  in  a  will  be  almost  symmetric  about  / 1  =0.5  (see  Figure  2).  The  coefficient  of  0  in  a  is 
the  product  of  / 2  and  a  term  symmetric  about  / 1  =  0.5  (see  Figure  3). 

If  Q2  =  0,  then  the  two  terms  in  (27)  involving  £3  are  replaced  by  a4«3  and  /1C3 
respectively,  and  the  term  involving  Qi  in  (28)  vanishes.  Equation  (30)  remains  as  is.  In 
(31),  the  two  terms  containing  Qi  in  the  square  brackets  are  to  be  replaced  by  2/2  and 
2/2/i  respectively.  The  symmetries  discussed  in  the  preceding  paragraph  still  hold. 

If  Q\  =  —j2ir2,  j  =non-zero  integer,  i.e.,  sinhQ2  =  0,  then  as  pointed  out  after 
equation  (A8),  C3  and  d3  are  infinite.  In  (C,  A),  d.$  is  multiplied  by  sinhQ2  so  ( C,A )  is 
bounded  when  the  basic  one-fluid  problem  is  at  neutral  stability,  making  the  numerator 
r  +  unbounded  due  to  the  C3~term.  In  this  case,  our  perturbation  procedure  is  not 

valid,  and  must  be  replaced  by  the  perturbation  procedure  for  a  double  zero  eigenvalue1. 
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Graph  of  a  versus  the  depth  of  the  lower  fluid.  axAT’  =  0.001.  R  =  1,  $  -  1,  a 
0.1,1.0,5.0,10.0,20.0. 
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Short  wave  asymptotics 

When  a  is  large,  the  boundary  conditions  at  z  =  0, 1  are  irrelevant  and  the  asymptotic 
behavior  of  o  is  the  same  as  in  the  case  of  solid  boundaries.  The  following  expansion  was 
derived  (Eq.(29).  Ref.  2)  by  scaling  the  z-variable  with  1/a  and  taking  the  distinguished 
limit  a  — »  oo,  So2  =  0(1): 


-( 

2a(±  +  l)[a,AT- 


,  .  1  ,  a25  \  1  . 

-  T  +0,?) 


'■s(s5F-'*‘,+»-7r)+0‘?-‘’>-  (32) 

If  Sa2  is  larger  than  0(1),  then  it  will  be  the  dominant  term  in  the  asymptotic  expansion 
but  the  other  terms  in  (32)  will  not  necessarily  be  the  correct  next-order  terms. 

The  numerical  calculations  checked  with  this  formula  for  / j  sufficiently  far  away  from 
0  and  1,  i.e.,  for  waves  that  are  short  enough  so  that  they  do  not  to  feel  the  presence  of 
the  boundaries.  The  coefficient  of  -  ^(JZa2)1/3  in  equation  (31)  behaves  as 


(1  -  *V 3)-(e~2l,Qo  -  e~2l3<is)  -h  e_2h<?3  _  e~2 
2 


-(1  +  *V3)^(e-2,1«>  -  e-2'^*)^  +  0(e"2a)  j 


as  a  — ►  oo 


so  is  exponentially  small  in  |aj. 


Long  wave  asymptotics 


For  long  wave  disturbances,  the  effect  of  surface  tension  is  0(a4)  and,  provided  Q 2  ^  0, 
a  ~  (Ba2/,  (  -  ^(15/J  +  30/2(/|  -  1)  +  7  +  <2(- 10  +  3/2)) 


'  6  +  h[f  +  *  '2  "  "*) 
+0(a4)  as  a  — ♦  0. 


•  v  *  •  •  •%*.••• 
•/  •  *  •  •  •  »  * 
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*  V  V  *.*  ’/’  *%  v  •*•»*•*.  •* 


When  Q%  —  0.  the  above  dominant  terms  are  0(ah)  and 

o  ~  to -  l\  -  1)5  -+  O(o6)  as  a  — *  0.  (35) 
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Thin-layer  effects 

In  the  linear  stability  analysis  of  parallel  shear  flows  composed  of  two  layers  of  fluids 
with  different  viscosities  and  densities,  the  effect  of  gravity  on  the  density  stratification 
can  be  countered  by  the  viscosity  stratification5.  A  linearly  stable  arrangement  is  possible 
with  the  more  dense  fluid  being  the  lower  fluid  if  it  is  also  very  much  less  viscous  than  the 
upper  fluid  and  if  that  layer  is  sufficiently  thin.  In  the  present  problem,  the  basic  flow  has 
no  shear,  but  we  find  that  the  stability  of  thin  layers  can  also  be  counter  to  intuition.  Here, 
the  role  of  viscosity  stratification  is  taken  over  by  the  thermal  conductivity  stratification. 
The  value  of  a  for  /j  small  is  essentially  =0)  where 

f>  -  o) 

Hence,  in  the  presence  of  a  thin  layer,  the  effect  of  surface  tension  and  the  differences  in 
density  and  cubical  expansions  are  dominated  by  the  difference  in  thermal  conductivity. 
The  coefficient  of  in  (36)  is  positive  for  Rayleigh  numbers  up  to  the  first  critical  value 
27tt 4  /  4  for  the  one-fluid  problem,  and  typically  looks  like  Figure  5.  At  the  first  critical 
value  of  the  Rayleigh  number,  the  coefficient  has  a  pole  and  is  infinite  at  the  critical  value 
of  a  =  ir /\'2.  As  the  Rayleigh  number  increases  above  the  first  criticality,  there  will  be  two 
values  of  o  for  which  the  one-fluid  problem  is  neutrally  stable3,  and  hence  the  coefficient 
will  have  two  poles,  until  the  Rayleigh  number  reaches  that  of  the  second  criticality,  when 
the  coefficient  will  have  three  poles,  and  so  on.  Our  perturbation  scheme  is  valid  away 
from  the  poles. 

For  small  /j.  a  <  0,  and  by  antisymmetry,  for  small  l?.  a  >  0.  Hence,  the  arrangement 
adth  a  thin  layer  of  a  less  dense  fluid  lying  below  the  more  dense  fluid  is  stable  to  long 
ind  order  1  wavelength  disturbances,  if  the  lower  fluid  has  the  lesser  thermal  conductivity. 


Real\(l  +  iV3)<2i  coth  Qj]  -  Q2coth  Qi  J  +  0(e2). 


'.v  -  '  v’.'  v  *^7 iw;.  ■.-:  '.-.  ■•.  'A1  -,v .r’AT w^t t  i*:- v.  «.-^r  >.-  «.■  n»  <r  ^  y» «-  g»  g».--i-? 

However,  for  short  wave  disturbances,  the  effect  of  the  thermal  conductivities  is  exponen¬ 
tially  small,  and  the  density  difference  and  surface  tension  dominate  the  stability  criterion. 
Figure  6  displays  o(a)  for  f  =  0.1,  <  =  100.0,  0  =  1.0.  S  =  10.0,  R  =  0.1  and  l\  =  0.05 
showing  linear  stability  even  though  Fluid  1  is  the  less  dense. 
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Graph  of  the  coefficient  of  -<<r  in  -  0)  versus  a  for  R  =  10.0,50.0,100.0,200.0. 

ie  amplitudes  decay  to  zero  for  large  a.  As  R  approaches  277r4/4,  the  peak  amplitude 
proaches  infinity  and  the  location  of  the  peak  approaches  a  =  7r/V 2. 
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Graph  of  a  versus  a.  f  =  O.i,  f  =  100.0,  0  ~  1.0,  S  —  10.0,  / j  =  0.05,  R  — 
0.1,1.0,10.0.  Fluid  1  is  the  less  dense  fluid.  At  R  —  0.1  and  1.0,  the  arrangement  is 
linearly  stable  for  all  a.  As  R  increases,  the  arrangement  will  become  less  stable.  For 
example,  at  R  =  10.0,  the  arrangement  is  unstable. 
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APPENDIX  :  Eigenfunct  ions  of  the  unperturbed  problem 


If  c=0,  the  variable  h  does  not  occur  in  the  right  hand  sides  of  (4),  (5)  ,(8),  or  in  the 
interface  conditions,  and  we  have  the  eigenfunction 


A  =  el 


The  adjoint  equations  yield: 


L0C  —  0, 


B'0C  =  0, 


where  C  =  (0' ,u‘ ,w' ,/T).  This  leads  to  the  equations1: 


£0’  +  RPw'  =  0, 


PAu^0--^-=O, 

dz 


du'  dw' 

~x~  +  "7TT  = 
ox  o  z 


Boundary  and  interface  conditions  are1: 


_  du' 

&'  =  w  =  — —  =0  at 
az 


at  2  =  /,,  f©  f  =  0. 


r^0’r  o 


ju*|  =  0. 


{juTf  =  0. 


,  du'  dw' 


1  dz  ^  dx^'  °’ 


•  *  «_*  •  •  -.*  s*  • 
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We  set  w'  —  u>(,e,az  etc.,  and  obtain  by  combining  the  equations: 


a2)V.  -  -  0. 


The  general  solution  of  this  equation  is 


W(,  =  c(  sinhQi^  +  c2  cosh  Q\z  +  c3sinh  Q2z  +  c4  cosh  Q2z 


+C5sinhQ3z  +  eccoshQ3z 


in  fluid  1,  and 


w'0  =  d\  sinh(Ji(z  -  1)  t  </2coshQi(z  -  1)  -  d3sinhQ2(z  —  1)  +  <f4CoshQ2(z  -  1) 


-d5sinh  Q3(z  -  1)  i-  d(jcosh  C?3(z  -  1) 


in  fluid  2,  where  Q\  is  the  complex  conjugate  of  Q3. 


Q]=  o2  +  (Ra2)l/\'*/3, 


Q7  =  a7  -  ( Ra7)'! 3, 


(?l=a2  +  (fio2),/3t-‘’/3. 


The  case  when  Q2  =  0  is  considered  later  (see  equations  A16-A18).  When  Q\  =  —  7r2, 


sinh(?2  =  0,  and  the  problem  reduces  to  that  considered  in  Ref.  1,  i.e.  the  perturbation 
of  a  double  zero  eigenvalue.  This  occurs  first  at  R  —  and  a  =  7r  \/2,  and  the  Q  i 
and  Q3  above  reduce  to  Q  j  and  Q2  of  equation  (B8)  of  Ref.  1.  In  fact,  sinh(J2  =  0  for 
R  =  ( j 2 7r 2  •+  a2)3 /'a2,  j  =non-zero  integer,  so  that  when  the  basic  one-fluid  problem  is  at 


neutral  stability,  our  c3  and  d3  are  oo. 


The  coefficients  ci  -  cc,  and  d\  -  d<_,  must  be  determined  such  that  the  boundary 


conditions  are  satisfied.  The  conditions  (A3)  at  the  walls  reduce  to 


=  5?“'“  =  J?w<’ =  °- 


•  .*•  .*•  v  »  .  v  •; 

V  s"  . 

^  m  1.*’.  *.*  *-•  *  • 


•  .  *  •  •  •  V  •  • 

•  •  •  •  • 
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At  z=0,  this  yields 


C2  +  c4  4"  Cg  —  0, 

Q]cz  +  Q\ca  4-  Q3C0  =  0«  (-^10) 

QU'2  ■+  Q$c4  +  qW  =  °- 

From  thi  *e  obtain  c2  =  c4  =  co  =  0.  At  z=l,  we  find 

<f2  +  ^4  4  df,  =  0, 


Qjd2  4-  Q\d^  4-  =  0,  (All) 

Q  jd2  +  Q2^4  4"  ^3^5  =  0. 


This  yields  d2  =  d4  =  <f6  =  0.  The  first  five  of  conditions  (A4)  lead,  after  eliminating  u' 
and  0*,  to  the  conditions 


ci  sinhQj/i  +  e3sinhQ2/i  4-  essinhQ3/i  =  -d\  sinhQt/2  -  d3sinh(?2/2  -  dssinh(?3/2, 


C1Q1  cosh  Q 1/1  4-  c3Q2cosh  Q2l\  +  c5Q3cosh  Q3/i  =  Q\d\  cosh(?j/2 
+  Q2^3  cosh  Q2/2  +  Qzd$  cosh  Q3/2, 

CjQjsinhQi/]  4- c3Q2 sinh Q2/i  4-  c5Q|  sinh  Q3l j  =  -Q\d  1  sinhQi/2 
-Q\d3  sinh  Q2/2  -  Q^efs  sinh  Q3/2, 
ciQ*  sinh  Q  j/i  4-  c3Q2sinh  Q2/j  4-  C5Q3  sinh  Q3/j  =  -Q\d\  sinh  Q]/2 

-  Q\d3  sinh  Q2/2  -  Q^d5  sinh  Q3l 2.  ( A13) 


>«**«•*  F.k as saLa^s^ar^ar^jg SJG3JCT 


r  TU’UJ^ 


ci(Qi  -  2o2Qi)coshQi/i  +  c3(Q|  -  2o2Q|) cashQ2^i  +  cs(Q|  -  2o2Q|)coshQ3/j 
=  d j (Qf  -  2o2Qi) cosh Qj/2  +  ^3(^2  ~  2a2Q2)  cosh  +  <^s(Q|  -  2a2Qa) cosh <33/2- 

From  these,  we  find  the  following  relations  that  will  be  useful  later: 

sinhQi/2 
Cl  'sinhQj/j' 


j  sinh  @2^2 

**  ~  ~d*  sinh  Q2I1’ 

sinh  Q3/2 
c  5  =  ~ds~r~r~p:  7  • 
sinh  Q3I1 


We  express  the  coefficients  in  terms  of  d$: 


ci  =  d[ 


C3  =  d> 


Q3  sinh  Q3  sinh  Q\l2  (l  +  »y/3) 
Q]  sinh  Q\ sinh  Qzl\  2 

Qasinh  Q3  sinh  Q2l2  (1  -  *v/3) 
Q2  sinh  Q2  sinh  Q3/1  2 

sinh  Q3/2 
C5~  5  sinh  Q3/1 


<fi  =  ~d$ 


<23sinhQ3sinh<2i/i  (1  +  *V3) 


Q\  sinh  Q\  sinh  Q3/1  2 

d  -  -d  Q3S*nhQ3s>nhQ2ii  (1  ~  «V5) 
3  5Q2sinhQ2sinhQ3/i  2 

When  a2  =  (iEa2)^3,  Q2  =  0.  The  general  solution  of  (A5)  is 


(*14) 


(A15) 


w'0  —  c  1  sinhQjz  +  c3z  +  cssinhQsz  in  fluid  1 


Wq  =  di  sinhQifs  -  1)  +  d$(z  -  1)  +  d5sinh<?3(z  -  1)  in  fluid  2. 


(*16) 


As  expected,  from  taking  the  limit  as  Q2  — *  0  in  (A  14), 

h  , 

C3  =  ~ra  3- 
*1 

Hence,  (A  15)  holds  for  cj,  C5  and  d\,  and 

Q3/2sinh(?3  (1  -  *V5) 

Cs~is  sinh Qa/"i  2~ 


(>417) 


(-418) 


.v« 

s’  •  ' 

vv  *. 
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